#include <stdlib.h>
#include <stdio.h>
#include <math.h>

#include "jpeglib.h"

static int
udiff(unsigned int const subtrahend, unsigned int const subtractor) {
    return subtrahend-subtractor;
}


static double
square(double const arg) {
    return(arg*arg);
}

/* Luminance, Chrominance macros. */

/* The following are weights of the red, green, and blue components
   respectively in the luminosity of a color
*/
#define PPM_LUMINR (0.2989)
#define PPM_LUMING (0.5866)
#define PPM_LUMINB (0.1145)

void
YCbCrtuple(int red, int grn, int blu,
               double * const YP, 
               double * const CbP, 
               double * const CrP) {
/*----------------------------------------------------------------------------
   Assuming that the tuple 'tuple' is of tupletype RGB, return the 
   Y/Cb/Cr representation of the color represented by the tuple.
-----------------------------------------------------------------------------*/
    
    *YP  = (+ PPM_LUMINR * red + PPM_LUMING * grn + PPM_LUMINB * blu);
    *CbP = (- 0.16874 * red - 0.33126 * grn + 0.50000 * blu);
    *CrP = (+ 0.50000 * red - 0.41869 * grn - 0.08131 * blu);
}


static void
psnr_color(int red1, int grn1, int blu1, int red2, int grn2, int blu2,
           double * const ySqDiffP, 
           double * const cbSqDiffP, double * const crSqDiffP) {

    double y1, y2, cb1, cb2, cr1, cr2;
    
    YCbCrtuple(red1, grn1, blu1, &y1, &cb1, &cr1);
    YCbCrtuple(red2, grn2, blu2, &y2, &cb2, &cr2);
    
    *ySqDiffP =  square(y1  - y2);
    *cbSqDiffP = square(cb1 - cb2);
    *crSqDiffP = square(cr1 - cr2);
}

static int
changed_yPsnr(int width1, int height1,
           double const ySumSqDiff, 
           double const crSumSqDiff, double const cbSumSqDiff,
           const char filespec1[], const char filespec2[]) {

    /* The PSNR is the mean of the sum of squares of the differences,
       normalized to the range 0..1
    */
    double const yPsnr = ySumSqDiff 
        / (width1 * height1) 
        / square(255.);

    double z = 10 * log10(1/yPsnr);
//    fprintf(stderr, "Y  color component: %f\n", z);
    if (z < 30.0)
        return 1;
     else
	return 0;
}

int frame_changed(char *filespec1, char *filespec2)
{
	struct jpeg_decompress_struct cinfo1, cinfo2;;
	struct jpeg_error_mgr jerr1, jerr2;;
	int row_stride1, row_stride2;
	JSAMPARRAY buffer1, buffer2;
	double ySumSqDiff, crSumSqDiff, cbSumSqDiff;
	unsigned int col;
	int retval;
	FILE *file1, *file2;
	unsigned char *p1, *p2;

	file1 = fopen(filespec1, "rb");
	file2 = fopen(filespec2, "rb"); 
	if (!file1 || !file2) {
		fprintf(stderr, "error opening files\n");
		return(1);
	}
	cinfo1.err = jpeg_std_error(&jerr1);
	jpeg_create_decompress(&cinfo1);

	jpeg_stdio_src(&cinfo1, file1);

	jpeg_read_header(&cinfo1, TRUE);

	(void) jpeg_start_decompress(&cinfo1);

	cinfo2.err = jpeg_std_error(&jerr2);
	jpeg_create_decompress(&cinfo2);

	jpeg_stdio_src(&cinfo2, file2);

	jpeg_read_header(&cinfo2, TRUE);

	(void) jpeg_start_decompress(&cinfo2);

	if (cinfo1.output_components != 3) {
		fprintf(stderr, "Not an RGB image\n");
		return(1);
	}
		
	if (cinfo1.output_components != cinfo2.output_components || cinfo1.output_width != cinfo2.output_width || cinfo1.output_height != cinfo2.output_height) {
		fprintf(stderr, "Images are different sizes/modes!\n");
		return(1);
	}
	row_stride1 = cinfo1.output_width * cinfo1.output_components;
//	fprintf(stderr, "output_width:%d output_components:%d\n", cinfo1.output_width, cinfo1.output_components);

	row_stride2 = cinfo2.output_width * cinfo2.output_components;
//	fprintf(stderr, "output_width:%d output_components:%d\n", cinfo2.output_width, cinfo2.output_components);

	buffer1 = (*cinfo1.mem->alloc_sarray)
		((j_common_ptr) &cinfo1, JPOOL_IMAGE, row_stride1, 1);
	buffer2 = (*cinfo2.mem->alloc_sarray)
		((j_common_ptr) &cinfo2, JPOOL_IMAGE, row_stride2, 1);

	ySumSqDiff = 0.0;
	cbSumSqDiff = 0.0;
	crSumSqDiff = 0.0;

	while (cinfo1.output_scanline < cinfo1.output_height) {
		(void) jpeg_read_scanlines(&cinfo1, buffer1, 1);
		(void) jpeg_read_scanlines(&cinfo2, buffer2, 1);
		p1 = buffer1[0];
		p2 = buffer2[0];
		for (col = 0; col < cinfo1.output_width; ++col) {
			double ySqDiff, cbSqDiff, crSqDiff;
			psnr_color(p1[col*3], p1[col*3+1], p1[col*3+2], 
			   p2[col*3], p2[col*3+1], p2[col*3+2], 
			   &ySqDiff, &cbSqDiff, &crSqDiff);
			ySumSqDiff += ySqDiff;
			crSumSqDiff += crSqDiff;
			cbSumSqDiff += cbSqDiff;
		}
	}

	retval = changed_yPsnr(cinfo1.output_width, cinfo1.output_height, ySumSqDiff, crSumSqDiff, cbSumSqDiff,
               filespec1, filespec2);

	(void) jpeg_finish_decompress(&cinfo1);
	jpeg_destroy_decompress(&cinfo1);

	(void) jpeg_finish_decompress(&cinfo2);
	jpeg_destroy_decompress(&cinfo2);

	fclose(file1);
	fclose(file2);

	return(retval);
}
